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Annual  Report  /  AFOSR  Grant  69-1808 

Seismic  data  will  alow  us  to  Infer  the  physical  properties 
of  the  interior  of  the  earth  only  if  three  conditions  are  met: 

a)  the  nature  of  the  release  of  seismic  energy  Is 
understood, 

b)  the  effect  of  propagation  through  an  arbitrary  trans¬ 
mitting  medium  can  be  calculated,  for  a  class  of 
media  spanning  the  possible  range  of  true  conditions, 
and 

c)  a  unique  relationship  can  be  found  between  the 
physical  properties  of  the  medium,  and  the  propaga¬ 
tion  of  stress  waves  through  that  medium. 

None  of  these  conditions  Is  satisfied  for  real  seismic 
sources  and  the  real  earth.  Therefore,  seismologists  Inter¬ 
pret  seismic  data  in  terms  of  simple  source  models  and  earth 
models  for  which  the  above  conditions  hold.  In  this  work  we 
have  extended  theclass  of  models  for  which  those  conditions 
hold  so  that  the  models  may  describe  reality  more  closely. 
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a)  Seismi c  Sources 

It  has  previously  been  assumed  that  an  earthquake  focus 
could  be  represented  by  an  edge  dislocation  along  a  plane 
fault  surface.  This  source  produces  displacements  identical 
to  a  double  couple  set  of  body  forces,  and  is  thus  referred 
to  as  the  double  couple  source  model.  The  double  couple 
source  may  be  described  in  terms  of  pure  shear  strain,  and 
involves  no  net  force  nor  net  moment.  Such  a  source  is 
consistent  with  the  notion  that  the  energy  released  in  an 
earthquake  is  identically  that  produced  by  relative  dis¬ 
placement  across  a  fault.  Because  displacement  is  known 
to  occur  across  fault  surfaces,  the  double  couple  model  is 
appealing  on  physical  grounds  for  surface  earthquakes. 

At  some  depth  in  the  earth,  almost  certainly  within 
the  crust,  the  shear  stress  needed  to  overcome  friction 
across  an  arbitrary  surface  exceeds  the  shear  strength  of 
the  rock.  Below  this  depth,  the  concept  of  a  "fault" 
becomes  meaningless,  and  the  double  couple  model  loses 
some  of  its  physical  appeal  for  those  earthquakes  occurring 
in  the  mantle.  A  competing  model  has  been  developed  (Knopoff 
and  Randall,  1970)  to  describe  the  effects  of  a  rapid  phase 
change.  Rapid  changes  in  density  and/or  bulk  modulus  pro¬ 
duce  an  isotropic  radiation  pattern,  with  radial  first  P- 
motions  and  no  first  S-motions.  A  rapid  change  in  shear 
modulus  would  result  in  a  double  couple  radiation  pattern  if 


the  initfal  stress  were  homogeneous  plane  shear;  however, 
given  an  axial  strain  field  before  the  event,  the  radiation 
pattern  would  be  that  of  a  linear  vector  dipole.  The 
change  in  axial  strain  resulting  from  the  phase  change 
requires  a  volume  change,  thus,  the  total  radiation  field 
would  be  the  sum  of  the  isotropic  radial  field  and  the 
linear  vector  dipole  field.  This  model  will  be  referred 
to  as  the  "Compensated  Linear  Vector  Dipole." 

In  Randall  and  Knopoff  (1970)  the  effects  of  propaga¬ 
tion  through  a  layered,  attenuating  earth  model  with  a 
crust  are  considered.  A  statistical  test  is  given  for 
discriminating  between  the  double  couple,  the  compensated 
linear  dipole,  and  the  simple  volume  change  models,  based 
on  amplitudes  of  P  wave  first  arrivals.  This  method  is 
applied  to  five  deep  focus  earthquakes.  For  most  of  the 
shocks,  the  double  couple  seems  to  be  the  dominant  source 
mechanism,  but  there  in  evidence  that  the  compensated  linear 
vector  dipole  plays  a  significant  role  (i.e.,  65%  double 
couple,  35%  CLVD). 

b)  Wave  Propagation 

We  have  made  progress  in  extending  the  complexity  and 
generality  of  earth  models  for  which  surface  wave  velocity 
dispersion  curves  can  be  computed.  Previously  these  earth 
models  have  been  limited  to  perfectly  elastic  layered  half 
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spaces.  We  have  significantly  Increased  program  efficiency, 
allowing  more  detailed  models  to  be  processed  in  a  given 
amount  of  computer  time,  we  have  generalized  the  theory 
so  that  attenuation  in  the  various  layers  may  be  considered 
explicitly,  and  we  have  developed  a  transformation  which 
allows  the  layered  half  space  programs  for  Love  wave 
dispersion  to  be  applied  to  layered  spheres. 

Each  point  on  a  dispersion  curve  represents  a  solution 
of  an  eigenvalue  problem  Involving  the  structural  parameters 
and  the  two  variables,  k,  the  wave  number,  and  uu,  the 
frequency.  The  several  methods  for  computing  dispersion 
curves  differ  only  in  their  representation  of  the  "dispersion 
function,"  a  matrix  product  which  enters  Into  the  eigenvalue 
equation  in  k  and  uu  .  Although  the  various  representations 
are  equal  when  fully  expanded,  the  speed  of  computation  and 
degree  of  physical  insight  provided  depend  strongly  on  the 
compressed  form  of  the  representation.  We  have  shown  that 
the  various  representations  fall  into  two  categories  (Schwab, 
1970) 

"The  Thomson-Haskel 1  technique  synthesizes,  or  builds 
up  the  surface-wave  dispersion  functions  by  constructing 
layer  matrices  which  relate  the  components  of  motion  at  one 
Interface  In  a  layered  structure  to  those  at  the  next,  ihe 
product  of  these  layer  matrices  then  relates  the  components 
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of  motion  at  the  deepest  Interface  to  those  at  the  free 
surface,  and  this  layer-matrix  product  Is  used  to  con¬ 
struct  the  dispersion  function.  Knopoff's  technique 
begins  with  the  Immediate  construction  of  the  dispersion 
function  In  Its  full  determlnantal  form,  and  then  analyzes, 
or  decomposes  the  determinant  Into  a  product  of  Interface 
matrices,  which  are  derived  from  submatrices  of  the  de¬ 
terminant.  Each  of  these  Interface  submatrices  relates 
the  components  of  motion  In  the  layer  on  one  side  of  the 
Interface  to  those  In  the  layer  on  the  other  side." 

"Randall  (1967)  performed  the  Initial  programming  of 
Knopoff's  technique  for  Rayleigh  waves  and  demonstrated 
the  validity  of  the  method.  During  the  recent  work  on 
the  optimization  of  the  computational  aspects  of  the 
Thomson-Haskel 1  formulation  (Schwab  and  Knopoff,  1970), 

It  became  evident  that  much  the  same  optimizing  processes 
could  be  applied  to  Knopoff's  technique,  and  that  they  would 
lead  to  considerable  simplification  In  the  computation  of 
the  Raylelgh-wave  dispersion  function,  as  well  as  Improving 
the  speed  of  computation.  Two  modifications  of  the  original 
form  of  Knopoff's  technique  bring  about  these  Improvements." 

"In  Its  Initial  form  the  method  employed  interface- 
matrix  elements  which  were  complex.  A  simple  change  In 
notation  allows  all  elements  to  be  expressed  as  either  pure 
real  or  pure  Imaginary  quantities.  Using  the  technique  given 
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by  Schwab  and  Knopoff,  this  allows  all  computations  to  be 
carried  out  with  pure  real  numbers.  This  eliminates  the 
extra  time  required  to  compute  matrix  products  when 
complex  elements  are  involved.  It  also  removes  the  problems 
inherent  in  bracketing  and  refining  roots  of  a  dispersion 
function  having  both  real  and  imaginary  parts,  since  this 
function  now  becomes  a  pure  real  quantity." 

"The  second  change  merely  requires  that  the  elements 
of  the  interface  matrices,  which  are  defined  as  fourth- 
order  determinants  for  Rayleigh-wave  computations  and 
second-order  determinants  for  Love  wave  computations,  be 
evaluated  analytically  before  numerical  computations 
begin.  In  the  final  representation  of  the  Rayleigh-wave 
dispersion  function  derived  from  Knopoff 's  method,  numerical 
evaluation  of  the  matrix  elements  directly  from  the  determl- 
nantal  definitions  leads  to  significant  loss  in  speed  of 
computation,  since  four  of  the  elements  in  each  interface 
matrix  are  seen  to  vanish  when  evaluated  analytically.  The 
expressions  for  the  other  elements  are  quite  simple,  thus 
considerable  computation  time  is  saved  by  elimination  of 
these  determinantal  definitions  from  the  actual  program." 

The  modifications  for  Rayleigh-wave  dispersion  compu¬ 
tations  were  shown  to  increase  the  speed  of  the  Thomson- 
Haskell  technique  by  a  factor  of  almost  60,  based  on  a 
test  in  which  the  original  and  modified  programs  were  used 
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to  compute  dispersion  curves  for  the  same  model  (17  layers 
plus  half  space)  on  the  same  machine.  This  improvement 
should  help  to  reduce  computation  costs  for  those  involved 
in  wave  propagation  studies  under  the  VELA  project.  The 
modified  Thomson-Haskel 1  method  and  the  Knopoff  method 
compared  as  follows  (Schwab,  1970): 

"Of  the  various  versions  of  the  Thomson-Haskel 1  for¬ 
mulation  for  Rayleigh  wave  dispersion  computations,  the 
reduced-delta-matrix  extension  Is  the  most  powerful,  i.e., 
the  fastest  which  contains  the  feature  controlling  the 
loss-of-precis ion  problem  occasionally  encountered  by  the 
original  formulation.  The  results  of  the  delta-matrix 
extensions  are  actually  contained  in  Knopoff 's  work,  which 
appeared  earlier  than  these  extensions,  and  these  results 
are  obtainable  directly  from  his  formulation  without 
recourse  to  delta-matrix  theory.  The  flexibility  of 
Knopoff 's  method  is  used  to  devise  a  new  representation  of 
the  Rayleigh  wave  dispersion  function  which  is  more 
powerful  than  the  most  powerful  of  the  Thomson-Haskel 1 
versions,  i.e.,  it  contains  the  loss-of-precision  control 
feature  and  is  about  38%  faster  than  the  reduced-delta- 
matrix  extension;  and  in  fact,  is  about  12%  faster  than  the 
fastest  of  the  Thomson-Haskel 1  versions." 

"Knopoff's  method  is  particularly  suitable  for  Rayleigh 
wave  dispersion  computations  involving  the  use  of  both 


solid  and  liquid  layers.  The  technique  for  the  construction 
of  solid-solid  interface  matrices  can  be  generalized,  in  an 
obvious  manner,  to  the  construction  of  sol  id- 1 iquid,  liquid- 
liquid,  and  liquid-solid  interface  matrices.  The  only 
difference  between  Knopoff's  method  and  the  Thomson-Haskel 1 
technique  for  the  computation  of  Love  wave  dispersion  is  in 
the  matrix  representations  of  the  dispersion  function.  The 
dispersion  functions  themselves,  except  for  a  possible  sign 
reversal,  the  characteristic  times,  and  the  actual  program 
speeds  are  all  identical. " 

The  introduction  of  complex  elastic  moduli  in  the  layers 
allows  the  effects  of  attenuation  to  be  considered  directly. 
Preliminary  results  show  that  for  a  frequency- independent  Q 
model  typical  of  those  in  the  literature,  a  reduction  in 
phase  velocity  of  less  than  10~^  km/sec  occurs.  This 
validates  the  usual  assumption  that  the  direct  effects  of 
attenuation  on  phase  velocity  are  not  detectable  seismically. 
The  method  also  allows  direct  computation  of  the  attenuation 
factors  for  Love  and  Rayleigh  waves  as  a  function  of  period. 
Preliminary  results  indicate  that  the  upper  mantle  acts  as 
a  "notch  filter;"  waves  of  certain  periods  are  almost  undamped, 
while  neighboring  periodsare  strongly  attenuated.  This  effect 
may  explain  many  of  the  puzzling  aspects  of  higher  mode  surface 
wave  observations,  such  as  the  propagation  of  Lg  waves  across 
continents  but  not  oceans. 
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Biswas  and  Knopoff  (1970)  have  developed  a  transformation 
which  allows  Love  wave  dispersion  curves  for  a  layered  half¬ 
space  to  be  converted  for  use  on  a  layered  sphere.  Although 
the  equations  of  motion  for  the  toroidal  oscillations  of  a 
sphere  may  be  solved  directly,  the  half-space  to  sphere 
conversion  offers  certain  advantages:  first,  dispersion 
curves  for  a  layered  half  space  can  be  computed  much  more 
rapidly,  by  the  methods  outlined  above,  than  those  for 
the  layered  sphere,  and  second,  the  conversion  method 
allows  calculation  of  the  phase  velocity  at  frequencies 
other  than  those  corresponding  to  the  Integral  toroidal 
modes.  Normal  mode  frequencies  have  been  calculated  from 
the  converted  phase  velocities,  and  are  indistinguishable 
from  those  obtained  by  the  exact  spherical  computations. 

The  ability  to  compute  dispersion  curves  as  a  continuous 
function  of  frequency  is  especially  Important  when  determining 
group  velocity  curves.  The  group  velocity  U  Is  obtained  by 
differentiating  a  function  of  the  phase  velocity,  c,  with 
respect  to  frequency: 

1/U  ■  d«u(w/c)/du )  . 

Where  c  Is  known  only  at  widely  separated,  discrete  points, 
differentiation  is  not  justified,  and  the  group  velocity 
will  not  be  known  even  at  those  discrete  points.  Interpo¬ 
lation  by  standard  procedures  can  help  only  slightly.  The 
earth  flattening  transformation,  however,  provides  a 
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continuous  phase  velocity  curve  from  which  the  group 
velocity  curve  may  be  obtained  by  differentiation. 

c )  Invers i on 

Theoretical  progress  in  the  inversion  of  non-linear 
systems  of  equations,  and  the  developments  of  rapid  dis¬ 
persion  curve  computations  mentioned  above,  now  permit 
us  to  find  the  range  of  reasonable  earth  models  which 
could  produce  seismically  observed  surface  wave  velocities. 
Without  introducing  arbitrary  restrictions  on  the  models 
considered,  unique  earth  structures  cannot  be  found.  How¬ 
ever,  the  range  of  possible  structures  consistent  with 
the  experimental  data  can  be  explored.  We  are  currently 
applying  the  so-called  "hedgehog"  technique.  A  criterion 
of  fit  is  selected,  consistent  with  the  accuracy  of  the 
data,  and  a  model  is  found  by  steepest  descent,  Monte 
Carlo  search,  or  guesswork  which  satisfies  the  criterion  of 
fit.  A  set  of  parameters  is  then  established  to  describe 
possible  models,  and  the  parameters  are  then  incremented  to 
produce  new  models.  Dispersion  curves  for  those  models  are 
then  compared  with  the  seismic  data,  and  the  process  is 
repeated  for  successful  models.  Thus,  the  "volume"  in 
parameter  space  of  the  satisfactory  models  is  systematically 
eaten  away,  as  implied  by  the  term  "hedgehog."  The  importance 
of  rapid  dispersion  calculations  is  evident. 
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We  are  presently  investigating  extensions  of  the 
inversion  scheme  which  will  determine  the  number  of 
equivalent  degrees  of  freedom  In  a  set  of  geophysical 
data  and  select  the  optimum  set  of  data  points  from 
continuous  data. 

d)  Appl icat ions 

The  surface  wave  dispersion  techniques  described 
above  have  been  applied  to  Rayleigh  waves  across  the 
East  Pacific  Rise.  In  each  case  the  "one  station 
method"  was  used,  requiring  that  the  earthquakes 
studied  have  vertical  strike  slip  sources.  The  intent 
of  these  studies  was  to  limit  the  propagation  path  to 
a  well-defined  region,  so  that  the  dispersion  curves 
would  not  represent  averages  over  quite  different  struc¬ 
tures.  For  this  reason,  rather  small  shocks  were  studied. 

The  earthquakes  were  in  the  Rivera  fracture  zone,  a  trans¬ 
form  fault  between  two  sections  of  the  East  Pacific  Rise, 
so  the  shocks  were  almost  certainly  vertical  strike  slips. 

The  stations  were  located  In  Culiacan  and  La  Paz.  Inversion 
was  performed  both  by  standard  least  squares  procedure 
(Knopoff,  Schwab  and  Schlue,  1970)  and  by  the  "hedgehog" 
method  (Knopoff,  Schlue  and  Schwab,  1970). 

The  path  from  the  Rivera  Fracture  Zone  to  Culiacan  was 
characterized  by  a  transitional  oceanic  structure  with  a  well 
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developed  low  velocity  zone  for  shear  waves.  The  lid 
velocity  was  around  4.4  to  4.45  km/sec,  while  the  channel 
velocity  was  41  km/sec.  The  top  of  the  low  velocity 
zone  occurred  at  about  50  km.  below  sea  level.  The  data  ’ 
did  not  extend  to  long  enough  periods  to  determine  the 
depth  of  the  bottom  of  the  low  velocity  zone.  The  Rivera 
Fracture  Zone  to  La  Paz  path  showed  rather  unusual 
character.  The  phase  velocity  decreased  slightly  with 
increasing  period,  in  the  range  25-65  seconds.  The 
least  squares  inversion  produced  a  structure  with  a  very 
thick  lid,  90  km,  and  an  extremely  low  channel  velocity 
of  3.5  km/sec.  Another  model  in  satisfactory  agreement 
with  the  data  has  a  decrease  in  velocity  with  depth  at 
50  km,  and  another  decrease  at  90-100  km.  The  velocity 
in  the  deep  channel  is  4.3  km/sec.  In  any  case,  there  is 
evidence  of  a  strong  anomaly  along  the  path  from  Rivera 
Fracture  Zone  to  La  Paz. 

Summary 

Significant  advances  have  been  made  in  the  degree  of 
sophistication  of  the  earth  models  and  fault  models  to  which 
standard  geophysical  analysis  may  be  applied.  An  alternate 
to  the  double  couple  fault  model  has  been  developed,  and  is 
called  the  compensated  linear  vector  dipole.  This  source 
model  is  appropriate  for  a  rapid  phase  change  in  a  region 
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of  uniaxial  stress.  A  statistical  test  has  been  devised 
to  distinguish  between  the  two  alternate  source  models. 
Several  deep  focus  earthquakes  have  been  examined  with 
this  test.  The  double  couple  source  appears  predominant, 
although  the  compensated  linear  vector  dipole  is  also 
significant. 

Tremendous  improvements  have  been  made  in  the  speed 
with  which  dispersion  curves  can  be  com  puted  for  layered 
half  space  earth  models.  Attenuation  in  the  layers  can 
now  be  included  explicitly.  A  transformation  from  the 
flat  layered  half  space  to  the  layered  sphere  makes  use 
of  the  fast  dispersion  calculations  for  the  half  space, 
and  allows  calculation  of  phase  velocities  at  frequencies 
other  than  those  of  the  normal  modes  of  the  earth. 

The  "hedgehog"  procedure  has  been  developed  and 
programmed  for  the  inversion  of  surface  wave  dispersion 
curves.  This  systematically  explores  the  set  of  earth 
structures  compatible  with  the  data. 

The  above  techniques  have  been  applied  to  data  for 
Rayleigh  waves  across  the  East  Pacific  Rise.  The  path 
between  the  Rivera  Fracture  Zone  and  Culiacan  appears 
similar  to  normal  oceanic  transitional  mantle.  The  path 
from  the  Rivera  Fracture  zone  to  La  Paz  seems  anomalous, 
with  an  ultra- low  velocity  zone  beginning  at  about  100  km. 
depth. 
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A  new  model  for  a  seismic  source  mechanism  has  been  developed.  The 
new  model  describes  the  action  of  a  rapid  phase  change  as  a  seismic 
source  and  is  an  alternate  hypothesis  to  the  double  couple  model. 

Great  Improvements  have  been  made  in  the  efficiency  of  programs  for 
calculating  surface  wave  dL&parsion  curves  for  a  layered  half  space. 
Attenuation  in  the  layers  may  be  taken  into  account  directly.  A 
transformation  has  been  developed  which  maps  Love  wave  dispersion 
curves  for  a  layered  half  space  into  layered  sphere  dispersion  curves, 
thus  allowing  use  of  the  more  efficient  half  space  dispersion  programs. 
Hedgehog  inversion  systematically  identifies  a  range  of  acceptable 
solutions  to  a  system  of  nonlinear  equations.  The  above  developments 
have  been  employed  to  analyze  Rayleigh  wave  dispersion  across  the 
East  Pacific  Rise, 
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